!The convergence test comes back before the place is properly setup..
!Need atleast Height+ width (3*Width) steps before convergence check
!Start managing with github.
!Parelise with openMP.

module inputs
implicit none

integer, parameter	:: dp=selected_real_kind(15,300)
integer, parameter	:: number_of_threads	= 2 
integer, parameter	:: Conv = 1e4
integer, parameter	:: Width = 100
integer, parameter	:: Max_N = 1e4
real(kind=dp), parameter	:: Thy_0 = 1000d0
end module inputs

Program Coax

!This program calculats the electrostatic potential in a wire
!The program should be able to use all the processors on a shared memory machene.
!Hopefully will be openMP 4.0 Compatable (gpgpu support)

use omp_lib
use inputs

Implicit none


integer	:: i,j		!array loops
integer	:: n 		!Itterations loop
integer	:: ierr 	!error checker
integer	:: Height	!Height of the grid
integer	:: Metal_Height	!Height of the metal
integer	:: Metal_Width	!Width of the metal


!Convergence variables

real(kind=dp)	:: Point_A, Point_B, Point_C	!The 3 points the Convergence is mesured at
real(kind=dp)	:: Point_A_Old=1d0, Point_B_Old=1d0, Point_C_Old=1d0	!The 3 points to check against
integer	:: Point_A_i, Point_B_i, Point_C_i	!The i Co-ord of the 3 points
integer	:: Point_A_j, Point_B_j, Point_C_j	!The j Co-ord of the 3 points


real(kind=dp), dimension(:,:), allocatable	:: grid 

open(unit=10,file="Results.log")
open(unit=11,file="Grid.dat")

!sets the number of threads
call OMP_set_num_threads(number_of_threads)

!Get the size of the width (Height is double this)
Height = 2 * Width
allocate(grid(0:width+1,0:height+1),stat=ierr)
if (ierr /= 0) stop "Error Alocating grid"
Metal_Width = nint(dble(Width)/10d0)
Metal_Height = nint(dble(Height)/4d0)



Point_A_i = 0
Point_A_j = nint(7.5d0 * 0.05d0 * dble(Height))

Point_B_i = nint(6.0d0 * 0.1d0  * dble(Width))
Point_B_j = nint(5.0d0 * 0.05d0 * dble(Height))

Point_C_i = nint(5.0d0 * 0.1d0 * dble(Width))
Point_C_j = nint(12.5d0 * 0.05d0 * dble(Height))


!set up initial conditions
call grid_setup(grid,Height,Width,Metal_Height,Metal_Width,Thy_0)

!Itteratation starts
itt:	 do n=1,Max_N	!Number of itterations

!$OMP PARALLEL

!Calculate the positions along the centre
i=0
!$OMP DO PRIVATE(j) 
Core:		do j=(Metal_height+1),Height
			grid(i,j) = Calculate_Core(grid(i,j+1),grid(i,j-1),grid(i+1,j))
		end do Core
!$OMP END DO

!$OMP DO PRIVATE(i,j)
!calculate the lower half of the metal
Bottom:		do j=1,Metal_Height
 		do i=Metal_Width,Width
			!Calculate the positions in the main non-metal
			grid(i,j) = Calculate_main(grid(i,j+1), grid(i,j-1), grid(i-1,j),grid(i+1,j)) 
		end do
		end do Bottom
!$OMP END DO

!$OMP DO PRIVATE(i,j)
!Calculate the top part of the metal
Top:		do j=(Metal_Height+1),Height
		do i=1,Width
			grid(i,j) = Calculate_main(grid(i,j+1), grid(i,j-1), grid(i-1,j),grid(i+1,j)) 
		end do
		end do Top
!$OMP END DO

!$OMP END PARALLEL

!Check for convergence
	if (mod(n,20)==0 .and. n>(2*Width**2)) then
	Point_A = grid(Point_A_i,Point_A_j)
	Point_B = grid(Point_B_i,Point_B_j)
	Point_C = grid(Point_C_i,Point_C_j)

	if (abs(Point_A_old-Point_A) < Conv) then
		if (abs(Point_B_old-Point_B) < Conv) then
			if (abs(Point_C_old-Point_C) < Conv) then
				write(10,*),"Converged at itteration ",n
				exit
			end if
		end if
	end if


	Point_A_old = Point_A
	Point_B_old = Point_B
	Point_C_old = Point_C
	end if


	end do itt
!Itteration ends


!Write results to file

rite:	do j=0,Height+1
	do i=0,Width+1
		Write(11,*),i,j,Grid(i,j)
	end do
	write(11,*), 
	end do rite



CONTAINS


SUBROUTINE grid_setup(grid,Height,Width,Metal_Height,Metal_Width,Thy_0)
!This subroutine fills the array with the initial conditions including the
!- bottom which is a logrithmic equation.

implicit none
integer	:: i,j
integer, intent(in)	:: Height
integer, intent(in)	:: Width
integer, intent(in)	:: Metal_Height
integer, intent(in)	:: Metal_Width
real(kind=dp), dimension(0:Width+1,0:Height+1), intent(inout)	:: grid
real(kind=dp), intent(in)	:: Thy_0



!fill the array with 0
grid = 0d0
!Calculate the bottom.
j=0
do i=(Metal_Width+1),Width
	grid(i,j) = Thy_0 * log(dble(Width)/dble(i)) / log(dble(Width)/dble(Metal_Width))
end do
write(10,*),'Bottom Calculated'


!calculate the inner metal.

do i=0,Metal_Width
do j=0,Metal_Height
	grid(i,j) = Thy_0
end do
end do

END SUBROUTINE grid_setup



FUNCTION Calculate_Core(up,down,right)
!This subroutine calculates the electrostatic potential along the core of the wire.
!Due to the radial symetry the left point = right point

implicit none
real(kind=dp)	:: Calculate_Core
real(kind=dp)	:: up, down, right

Calculate_Core = 0.25 * up + 0.25 * down + 0.5 * right 

END FUNCTION Calculate_Core



FUNCTION Calculate_Main(up,down,left,right)
!This calcualtes the electrostatic potential for the majority of the wire.
!It does this by finding the average of the surroundings

implicit none
real(kind=dp)	:: Calculate_Main
real(kind=dp)	:: up, down, left, right

Calculate_Main = 0.25 * up + 0.25 * down + 0.25 * left + 0.25 * right

END FUNCTION Calculate_Main

END PROGRAM Coax
